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PREFACE 


The  HEMP  3D  numerical  technique  for  calculating  elastic-plastic  flow  in  three 
space  dimensions  and  time  was  developed  in  the  late  1960s  at  the  Lawrence  Livermore 
National  Laboratory  with  funding  from  the  Material  Science's  Office  of  the  Defense 
Advance  Research  Projects  Agency  (DARPA).  The  numerical  scheme  programmed  for 
the  CDC  6600  computer  by  John  French  was  first  presented  at  the  Second  International 
Conference  on  Numerical  Methods  in  Fluid  Dynamics  at  the  University  of  California, 
Berkeley,  September  15-19, 1970.  The  usefulness  of  the  computer  simulation  program  at 
that  time  was  limited  by  the  lack  of  adequate  three  dimensional  graphics.  With  continued 
funding  from  DARPA,  a  production  program  was  developed  for  the  CDCSTAR 
Computer  by  Eugene  Cronshagen  using  vector  programming,  including  three 
dimensional  graphics.  This  work  was  published  in  1975  as  UCRL-57574,  "A  Method  for 
Computer  Simulation  of  Problems  in  Solid  Mechanics  and  Gas  Dynamics  in  Three 
Dimensions  and  Time."  The  program  has  for  many  years  been  operating  on  the  CRAY-1 
computer  at  the  Lawrence  Livermore  National  Laboratory.  Presented  here  is  an  update  of 
the  1975  report  that  includes  the  sliding  surface  routines  programmed  by  Robert 
GuUiford. 
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HEMP  3D — A  Finite  Difference  Program 
for  Calculating  Elastic-Plastic  Flow 

INTRODUCTION 


The  HEMP  3D  program  can  be  used  to  solve  problems  in  solid  mechanics 
involving  dynamic  plasticity  and  time  dependent  material  behavior  and  problems  in  gas 
dynamics. 

The  equations  of  motion,  the  conservation  equations,  and  the  constitutive  relations 
listed  below  are  solved  by  finite  difference  methods  following  the  format  of  the  HEMP 
computer  simulation  program  formulated  in  two  space  dimensions  and  time.1 

A.  Equations  of  motion. 


0 

dx 

P —  = 

+^k 

Hdt 

dx 

dy 

dz 

iQ 

p*z- 

dT^ 

dTn 
+  -Z- 

H  dt 

=  dx  ' 

dy 

dz 

iii) 

dz 

P —  = 

dr  dr„ 

+-2L-1 

Hdt 

dx 

dy 

dz 

B.  Conservation  of  mass  and  energy. 

i)  ^  =o  ;  M  =  Mass  element 

dt 

C.  First  law  of  thermodynamics. 

i)  E  =  -(P  +  q)V  +  V[sja  +  s„e„  +  saea  +  +  T +  Taen] 

E  =  Internal  energy  per  original  volume 
V  =  Relative  volume  =  ptyp 
p  =  actual  density 
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pO  =  reference  density  of  equation  of  state 


D.  Velocity  strains. 


i) 


ii) 

iii) 


E.  Stress  deviator  tensor. 


iv) 

v) 


vi) 


dy 

dz  +  dy) 


F. 


|i  =  shear  modulus 
Pressure  equation  of  state. 

i)  .  P  =  a(rj-l)  +  b(T]-l)2  +c(tj-1)3  +drjE 

ii)  tj  =  —  =  p  /  p°,  where  a,  b,  c  and  d  are  equation-of-state  constants 


G.  Total  stresses. 

i)  'LIX=-(f3+q)+s~ 

ii)  l„=-(P  +  q)  +  s„ 

iii)  =-(^ +  ?)  +  *„ 
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H.  Artificial  viscosity. 

i)  q  =  dpL2^j +  CLpLa^  . 

Co  and  Cl  arc  constants 

ii)  —  =  rate  of  strain  in  the  direction  of  acceleration 
dt 

L  =  measure  of  grid  size 
a  =  local  sound  speed 
p  =  local  density 

I.  von  Mises  Yield  Condition 


Y  =  plastic  flow  stress 
y =a{b+e')’ 

ep  =  equivalent  plastic  strain 

21  =  second  invariant  of  the  deviatoric  stress  tensor 

a,  b  and  c  are  flow  stress  constants. 


3iuiv  siqiaa  a  a  •  w.sk»  li  *  i  as  i 


The  finite  difference  equations  that  integrate  the  equations  of  physics  are  based  on 
the  divergence  theorem.  Partial  derivatives  are  evaluated  by  summing  the  flux  through 
the  surface  enclosing  an  element  of  mass. 


Thus,  both  the  physics  equations  and  the  method  of  solving  them  are  in  a  conservation  form. 

The  physical  object  is  divided  into  zones  defined  by  eight  grid  points,  Rg.  1.  The 
grid,  (i.jjc)  moves  with  the  material  and  the  mass  within  a  zone  remains  constant.  In  the 
notations  that  follow  a  superscript  refers  to  the  time  centering  or  a  parameter  or  equation 
and  the  subscript  refers  to  the  space  centering. 


Three  vectors  are  associated  with  each  of  the  eight  grid  points,  g,  shown  in  Rg.  1. 
2  =  I 


f-4-r 


- 


Vector 

Components 

\ 

-A 

i 

\! 

-li 

4: 

(a;),  =  4  -^K),  =  *4  -  V 

i  i  • 

r  A  U  '  1 

44--^ 

fo),  =x2-x;(by)i  =y2-y1;(bi),  =  z2-21. 

1  r  2 

C: 

te).  =Jts-4/),  =y,-y»;(^X  -*»-*!• 

Vector 

Components 

X: 

(a;)2  =  *,  --t,;(a;)2  =y,  -y2;(a*)2  =  z,  -z,. 

r  *  c  i 

1  • 

!  >-”-^3 

B: 

fo)2 =*3-  -v(*/)2  =  y~\ -yAbt)7  =  h~h- 

X  2  8 

C : 

(c.)2  =  r»  -  Mc/),  55  y6-yAck)2  =  z6  -  V 

-1“ — >T7 

Veaor 

Components 

r— • + — <  ^ 

|  IB*'  • 

1  4-bi3 

X: 

(aA  =  *=  -  ^(a/)3 = y2-y3;M3  = 

B: 

(b,)3  =x4  -.r,;(b;)3  =  y4-y3;(b,)3  =  z4-zr 

IcT - A 

2 

C : 

(c-  \  =  xr~  4’(cy)3  =  y2  ~y3*{ck)3  =  Zy  -  z3. 

O 

Vector 

Components 

*  K),  =  ^ )4  =  * -y*'»iak)4  =  z,  -  *4- 
&  =  T.  -.W,(bj)4  =  y,  -y4;(b4)4  =  z,  -  z4. 

C-  k l  =  -t,  -  .ct;(cy)  =  y,  -  y4;(c*)4  =  zt-z4. 
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-► 

C 

2  =  6 


Vector 


Components 


i.i 

M 

- J 

-4:  (*X=x*-xAa/)s 

!  A !  ! 

I  ^ 

B:  (bi)1=xs-xs;(bJ)s 

L' _ u" 

i 

C:  fa)s  =  x:  ~xs'(cj)s 

Vector  Components 

*  (<04 = - ^;(a;)6 =>7- y*;(a*)4  =  *t-v 
B-  {b)6  =  *5  - *6;fe)6 = >s -y»(K\  =  *s  -*«• 
C:  fo)6  =  *2  - -t6;(c;)6  =  y2- y6\{ck)6  =  z2- z6. 


g«8 


Vector  Components 

*  (al)7  =  xt-xr;(ay)7*y,-y7;(al)7  =  2l-z7. 

fl:  fo)7  =*«- = y6  - y7;(^)7  ■*«-*?• 
C-  (c.)7 =*3- x,;(cj )7 =y3~  y7;(c*)7 -  -  v 


—  8 


B 


— r 

i 

*c  ' 


U--- _ I*' 


i-~J 


Vector  Components 

A:  (a,),  =  Xj  -  *,;(«/),  =ys -*;(<»*),  =  r5  -  z,. 

fl:  fa), = *7-*,;fa,), = y7 -y«;fa*), ^h-h- 

C-  fa),  =  *« -*;fa),  =  y<-yti{ct)t  =  z4~V 


f 
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Refer  to  Fig  1. 


Calculation  of  the  Volume  of  Zone  ©.  v@ 


Refer  to  Fig.  1. 

v© 


f BXA-C]  n  = 

L  JP  =  1 


bt  bj  K 

\ai  <*j 

C,.  Cj  ck 


=  [bi(ajCk  -  atc,)  -  bt{afik  -  akq)  +  bk(aicj  -  a/:,)]' 


i-i 


Repeat  for  g  =  2-4  8 


Calculation  of  the  Mass  of  Zone  1. 


p°  =  reference  density 
V°  =  initial  relative  volume 

v°  =  actual  volume  calculated  from  the  coordinates  at  time  t  =  0. 


8 


f  nQ\ 

=  j^£_  where  is  the  volume  at  time  t  =  n  and  Vq  is  the  relative 


volume.  Similarly, 


■ta 


ImV® 


v^t1  where  the  volume  vn+1  is  calculated  from  the  coordinates  at 


time  n+1. 


V01/2  =  —  (V0n+1+  V0n)  definition  of  relative  volume  at  t  =  n  +  1/2. 


The  following  acceleration  equations  are  applied  to  point  0  in  Fig.  2. 


Wi  j  Jc  =  -  [M©  +  +M@  +M@  +M©  +M©  +M©  +M@] 


f&\  1  \ f  arv  t 

ydthjjt  P"mL  *  % 


where 


1  ggJ  -  1 

P  A  Aj.k  40..m 


{  (^xx)©  [(yyj"yy)(zjY"zy)_(zyj*zy)G,|v'yv)] 
(^xx)©  [(yjj"yy)(zyj.zy)_(z]]"zy)(yyj'yy)] 

+  (^xx)©  [(yiv-ym)(zvrziii)-(ziv-ziii)(yvi-ym) 
+  (^xx)©  [(yvi-yin)(zii'zin)"(zvi-zin)(yn"yni)l 

"*■  (^xx)(§)  [(yy'yj)(zjy_zj)-(zy-zj)(yjy_yj)] 


f 
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+  (£xx)@  Kyn-yjKzv^^^n^^^v'yi^ 

+  Gxx)®  [(yi-ym)(ziv'ziii)-(zrzin)(yiv‘yii^ 

+  Oxx)@  [(yn-yni)(=i-ziii)-(zn‘zin)(yryin)]  \<jtk 


Fo  form  (1/p  9Txy/9y •)"*  ,  replace  each  Xxx  in  the  above  expression  with  TXy,  every  y 
with  the  corresponding  x  and  every  z  with  the  corresponding  y. 


The  x-direction  velocity  at  n  +  1/2  and  positions  at  times  n  +  1  and  n  +  1/2  are: 


•  *+1/2  _  - *-1/2  ,( dx}  A,* 


*+1/2 


_.*+l/2  _  1  /  «+l  ,  _«  \ 


Motion  in  tfig  y  Pirectifla 

dyY  1  r^T„  ar.T 

^J,r^[^+^+^L;wherc 


'iar- 


VP  *  J,y  4 
by  the  corresponding  value  of  Txy. 


=  same  as  (1/p  dXxx/3x)*4  ,  defined  above,  except  replace  each  Xxx 


f  ^  ^  ^  \n 

21 1  =same  as  (1/p  3Txy/3y^  t  ,  defined  above,  except  replace  each  Txy 


vp  ki .* 

by  the  corresponding  value  of  Xyy 


—  ■  =same  as  (1/p  cdzx/dz)]^  ,  defined  above,  except  replace  each  Tzx 


VP  *  kj.A 
by  th#*  corresponding  value  of  Tyz. 
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The  y-direction  velocity  at  time  n  +  1  and  positions  at  times  n  +  1  and  n  =  1/2  are: 


•  »+i/2  _  .-,0-1/2  ( dy\  „ 

y.j.t  - yij*  +!  ~r  I  & 
V  dt  Jijjt 


v"+1  =v"  4- v"+l/2Af*+,/2 


Motion  in  the  z  Direction 


ydt)iJik  p*jk 


a*. 

L  <?y  <9z 


I. 

-J*.y.* 


i  ar) \ 


VP  d*  Ay.* 

Zxx  by  the  corresponding  value  of 


-same  as  (1/p  3Xxx/9*)  ,  defined  above,  except  replace  each 


1  dT  X 

1  -same  as  (1/p  8Txy/3y)  ,jt  ,  defined  above;  except  replace  each  T 


VP  dy  Ay.* 

by  the  corresponding  value  of  Tyz. 


*y 


( 


1  9Z  Y 

—  =same  as  (1/p  ,•  ,  defined  above;  except  replace  each 


*  Ay.* 

Tzx  by  the  corresponding  value  of  Zzz- 


The  z-direction  velocity  at  time  n  +  1/2  and  positions  at  times  n  +  1  and  n  +  1/2 
are:  % 


^+1/2  -0-1/2  f  dzY  -o 

•-/-*  z«.y.*  +  j,  m 

V  ot  /  i,y.* 


o+l  _  o  ,  -0+1/2  A  ,0+1/2 
Zi.j,k  ~  Zi,j.k  +  Zi,j,k 


f 


11 


+  z, 


)• 


Calculation  of  Incremental  Strains 

The  finite  difference  mapping  procedure  to  calculate  the  surface  integral  of  zone 
©,  Fig.  1,  covers  the  surface  in  units  of  triangles.  The  velocity  associated  with  a  given 
triangle  is  taken  as  the  average  of  the  velocities  defined  at  the  triangle  comers.  The 
triangular  surface  area  vectors  are  calculated  to  point  out  of  the  zone  surface.  The  dot 
product  of  the  area  vector  with  the  direction  vector  multiplied  by  the  average  velocity 
gives  the  velocity  flux  through  the  surface  in  the  given  direction.  The  mapping  procedure 
actually  covers  the  zone  surface  area,  Fig.  1,  two  times.  The  difference  equations  used  to 
calculate 

~  ,  ^  and  ^  arc  given  explicitly  below. 

dxdydz 


The  remaining  velocity  derivatives  required  to  calculate  the  components  of  strain  arc 
calculated  by  replacing  x  in  those  equations  by  y  and  then  by  z  so  as  to  complete  the  set: 


dx 

dx 

dx 

dx 

dy 

dz 

dy 

dy 

dy 

dx 

dy 

dz 

dz 

dz 

dz 

dx 

dy 

dz 

Velocity  Derivatives  Corresponding  to  Zone  ©.  Fig.  1 . 


12vJ1/2 


where 
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(*A*)f=1  =  (*i  +X2+  xA),{xCA)t=l  =  (i,  +  i4  +  x5),(iflC)<=i  =  (i,  +  ij  +  x5). 


1  0  0 


fi  ai  ak 
b,  bj  bk 


=[(aA- "A  )],,,• 


*=i 


( Cxi‘\ 


1  0  0 

C:  C.  Cl 


\a:  aj  ak 


=  [(cA-cia.)]|=1. 


ii=i 


1  0  0 

(BXC  i)i  t  =  b,  b,  bk 

c,  ct  ct 


=((*A-V,  )],„■ 


I,=l 


The  above  steps,  written  for  g  =  1,  must  be  repeated  for  g  =  2  — >  8. 


dx 


■+I/2  f  \ 


£  •  ] +ia(CXA) ■  ] + 4(»c)- 


where 


and  are  as  defined  above. 


iAxSA.r 


0  1  0 

a,  a-,  ak 

b,  b ,  bk 


=KaA 


I*=i 


[£x'AA.r 


0  1  0 


C(.  Cj  ck 
I a:  dj  dk 


=  [-k«*-cta,)] 


f*l 


1=1 
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|0  1  0 

(ix£  A.rh  b-  fc* 


Pi  Cj  ck 


=  [-(Vi-V,)]f=1- 


l,=l 


The  above  steps,  written  for  g  =  1,  must  be  repeated  for  g  =  2  — >  8. 


^V+1/2 


Jy)® 


r 


i 


i12v®  j 


\  g 

r  £  [M*  X B)  j  +XC,(CXA)  J  +  itc(BXC)  ]]"'n 


where 


30(1  (*«c)t=t  arc  defined  above. 


(izi-f)  - 


0  0  1 


P.  «y  A* 


b,  bj  bk 


=[(«.*/ -<’A)],.1- 


1**1 


pxi.f)  - 


0  0  1 


\Ci  Cj  ck 


I A.  A;  At 


>1=1 


0  0  1 

h  bJ  b* 


C:  Cj  Ck 


=[(t,.c, -¥.)],„■ 


1=1 


The  above  steps,  written  for  g  =  1,  must  be  repeated  for  g  =  2  — »  8. 
o>y 

-r^-=  same  as  dx/dx  except  replace  x  by  the  corresponding  y 


— =  same  as  dx/dx  except  replace  x  by  the  corresponding  z 
dx 
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_ 


dy 


=  same  as  dx/dy  except  replace  x  by  the  corresponding  y 


dz 

— =  same  as  dx/dy  except  replace  x  by  the  corresponding  z 
9y 


dy 

-r-=  same  a s  dx/dz  except  replace  x  by  the  corresponding  y 
dz 


dz 

—  =  same  as  dx/dz  except  replace  x  by  the  corresponding  z 
dz 


Incremental  Strains 


. 


/ \»+l/2 


K>; 


>.+1/2 


'<D 


\dzj®  {dy)® 


( AVV+’/2  /  \*+l/2  /.  \«+r/2  .  yi+l/2 

[ir L  ~(A£»)(D  +  (Ae»)<j)  +(A£»)<i> 
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m?’ ■  w&> jsn  -  lif  ■£"} f 
(°x - 

(^>®  =w®  ^[(^r 4(v)r]+(a:)®' 

(r«)®=(r.)®+^.)®'1+(^)®' 

Note:  The  terms  8  that  have  been  added  to  the  stress  deviators  are  corrections  for  zone 
rotations.  If  a  zone  has  rotated  during  the  time  interval  Atn+1/2  =  tn+1  -  tn  the  stresses  at 
time  tn  must  be  recalculated  so  that  they  will  be  referred  to  the  x,  y,  z  coordinate  systems 
in  their  new  positions.1  The  correction  terms  are  given  by 

8nxx=-2(0HtT*+2(o]T*, 

=  +2  0)"  7^ -2  o>"  T£ , 
s:=+2m;T;-2  a;  r; = s-„  -  , 

K,  =  (C  -*"„)+ a"  K  ~  <*>."  C 
s; = ©;  (s; -s"„)+ ©;  r;  -  »;  r;, 
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5"  =0)*  fs"  -s"  )  +  g)',7’*-0)"T" 

zx  wy  y’z z  xx /  '  x  xy  z  yz 


where 


0)  =- 
*  2 


dy  dz) 


Pressure  Equation  of  State 


where  A  and  B  are  functions  of  the  volume  V  and  E  is  the  internal  energy. 


Total  Stresses 


von  Mises  Yield  Condition 

2'®  =  {[(O2 + (vf + (O2] + 2[(r„)2 + (tJ + (r„)2]}”'  • 
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t'rt+l  o  rrt+l  ^  ( v* 

*(D  “"a*  ""  3I  ®;  • 

If:  ATq1  <  0  use  the  stress  deviators  as  defined  above. 

If:  AT,**1  >-  0,  then  multiply  each  of  the  stresses  (s^)^ ,  (vWW  «d 


s,  =  2«J-a/3cos<p/3, 
s  2  =  2*j—a/3  cos(<p/ 3  +  2ji/ 3), 
s3  =  2*J-a/3cos(<t>/3  +  4^/3); 


where 

a=[(v«  +s»s» +sx*st,)-(Ti+TZ+TZ% 

b  =  V-  +  s~Ti  +  S°TZ  +  ~  2Tyzr*yT**\ 

/ 

s,  =  max  s,  ,s2>s3 
V  ~  ~  ~ 

/ 

s3  =  min  s,  ,s2,s3 
<  ~  ~  ~ 

s2  =  —  (s,  +  s3) 
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Artificial  Viscosity 

An  artificial  viscosity  is  required  to  permit  shocks  to  form  in  the  grid.  The 
artificial  viscosity,  q,  used  here  is  composed  of  a  quadratic  and  linear  function  of  the  rate 
of  strain.  The  quadratic  portion  is  a  generalization  to  three  dimensions  of  the  one 
dimensional  von  Neumann  q  for  calculating  shocks.2  The  linear  portion  provides 
damping  for  oscillations  that  can  occur  behind  the  shock  with  the  q  method  of  calculating 
the  shock  front  The  term  ds/dt  used  in  the  q  calculations  here  is  the  rate  of  strain  in  the 
direction  of  acceleration.3 


Ax,  Ay  and  Az  are  the  x,  y,  z  components  of  acceleration  respectively. 

L  =  measure  of  the  zone  size  taken  here  as:  3] zone  volume 
a  =  JFTp 

Co=  ~  2 
CL=  =1 


The  q  is  added  to  the  pressure  P. 
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Tensor  Artificial  Viscosity  for  Stabilizing  the  Grid 

For  quasi-static  problems  in  solid  mechanics  nonphysical  numerical  oscillations 
can  occur  in  the  grid  under  certain  boundary  conditions.  A  tensor  viscosity  based  on  the 
rate  of  strain  of  volume  elements  formed  by  the  zone  comers  is  used  to  damp  this  type 
oscillation.  Referring  to  Fig.  2  it  is  seen  that  surrounding  point  0  there  are  eight 
tetrahedrons  defined  by  the  comers  of  the  eight  zones.  A  Navier-Stoke  type  tensor 
viscosity  based  on  the  rates  of  strain  of  the  tetrahedron  volumes  is  calculated  for  each 
tetrahedron  that  contains  0,  Fig.  2.  The  details  for  calculating  the  components  of 
viscosity  for  the  tetrahedron  in  zone  ©  are  given  below. 

The  tetrahedron  corresponding  to  zone  ©  is  shown  in  Fig.  3.  The  grid  numbering 
follows  the  scheme  shown  in  Fig.  1.  Here  grid  point  1  corresponds  to  point  0  of  Fig.  2. 
The  finite  difference  integration  mapping  procedure  is  applied  to  the  four  surfaces  of  the 
tetrahedron  formed  by  vectors,  A,B,C,  of  Fig.  3. 

Volume  Vabc  formed  by  the  vectors  A,B,C,  of  Fig.  3  is: 

The  notation  for  the  components  of  the  vectors  is  the  same  as  used  for  the  vectors  of  the 
volume  of  zone  © . 


Velocity  Derivatives 


Velocity  derivatives  corresponding  to  the  tetrahedron.  Fig.  3. 


( Y,+1/2  (  ,  V 

l*  J  =  1 JM*  XB)i  +  xa(cXA)J  +  xK(BXC)T  +  xCD(EXD)i 


where 
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=(x,  +i2+i4);  xCA  =  (i,  +x4  +i5);  iBC  =  (i,  +x2  +i5);  x^  =  (i2  +i4  +x5); 


and  v~‘/2  =  +  v^J. 


This  expression  can  be  simplified  by  expressing  vectors  D  and  E  in  terms  of  vectors  A 
and  B  . 


where 


| AX  B-i)  =  [cijbk  - akbj );  (CXA-i)s ~(cjak  - cka^, and (flXC-i)  =  ~if>jck ~ bkc,\ 


n«+l/2 


X  -x,)(cXA)J+(x, -x4){BXC)-irm, 

ibc  r 


where 


AX B  j  =  -(a,b*  - afy);  (CXA-J)  =  -{cflk  - cka), and(BXC  j)  =  -{bfik  “  Vi> 


v.n«+i/2  ( 


X  §)■  k  +(i,  -i,)(C  X  A)  k+(x,-x,p  X 


where 


A  XBk  =  («,<>, -afy)-,  CXA-k-fa-cfl);  B  X  Ck  =(b.ct -bfi). 

ft  ft 

~  and  —  are  calculated  in  the  same  way  as  dx/dx,  but  replace  x  by  y  and  then  z. 

dx  dx 

ft  ft 

77-  and  —  are  calculated  in  the  same  way  as  dx/dy,  but  replace  x  by  y  and  then  z 

ay  ay 
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v 

—  and  —  are  calculated  in  the  same  way  as  dx/dz,  but  replace  x  by  y  and  then  z  . 
dz  dz 


Rates  of  Strain 

Components  of  the  rate  of  strain  of  the  tetrahedron  defined  by  vectors  A,  B,  C,  Fig.  3  are 


v  _  dx  |  cjy  t  dz 
v  dx  dy  dz 


Artificial  Viscosity 

Tensor  artificial  viscosity  for  tetrahedron  A,  B,  C,  Fig.  3,  is 


o+l/2  »+l/2  o+l/2 

?"‘,2  =  2ft[i=-“];  <"  =  2ft[i„-~];  C'2  =  2ft[£„-~]; 


o+l/2 


o+l/2 


o+l/2 


<1/2 = ^  =  a[4]  • 


where  pi  = 

c  (£ 

.  )  \jVAI)C 

L  ^ 

w©  J 

Cns  =  constant  =  10'2, 

pO  =  reference  density  of  zone  ©, 

V  =  relative  volume  of  zone  ©, 

The  above  components  of  the  tensor  artificial  viscosity  are  added  to  the  corresponding  components  ol 
stress  tensor  defined  at  time  n  +  1. 
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Increment  of  energy  dissipated  by  the  tensor  artificial  viscosity 
8 

AW©  =  +  qn£n  +  qja  +  q^  +  qja  +  qjy,].  At 

i- 1 

Here  i  =  1— >  8  are  the  eight  nodes  that  define  zone  ©. 


Energy  Equation 


Distortion  Energy  Increment 

Az£,/2  =  Ae  +s  Ae  +s  Ae  +T  Ae  +T  Ae  +T  Ae  1!!‘/2 

®  ®  L  «  «  yy  yy  «  xx  *  zy^^xy  Ayx^^yx  * 

sxx,  etc.  and  Ae^,  etc.  are  the  components  of  the  stress  tensor  and  increments  of  strain 
respectively  defined  at  the  zone  center. 


Total  Internal  Energy  per  Original  Volume 


r»+l  _ 

fc©  “ 


(  E*  -  |^[a(V"+1  )  +  />"]  +  •  (V-  -  V")  +  A z"+I/2 


l  +  -[fl(V"+1)](V"+,-V") 


AD 


i=\(rx,1+rm) 


Note:  It  has  been  assumed  here  that  the  pressure  equation-of-state  has  the  form  P  =  A(V) 
+  B(V)E. 
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Plastic  Strain 


In  the  following  definitions  of  plastic  strain  the  stress  deviators  at  time  (n+1)  are 
taken  as  the  values  after  the  yield  condition  has  been  satisfied.  If  yielding  has  not 
occurred,  these  equations  are  bypassed. 


Components  of  Plastic  Strain  Rate 


n  •  —  o»+l/2 

P  £<1+1/2—  £„ 


Pg«+I/2  _  £*+1/2  . 


P£*+l/2  _  £»+I/2 


p  £*+1/2  _  £*+1/2  . 


1 

r  cis--s„ 

Af«+ 1/2 

hi 

1 

re-*-*. 

Ar«+ 1/2 

2/i 

1 

>♦»  _  c"  -X 

Ar"+,/2 

L 

1 

xy  xy  xy 

Ar"+  ■i,2 

M 

1 

T^1  -T* 

A/*+ 1/2 

M 

1 

/y,  iy,  dyj 

Ar*+1/2 

/*+!  i/* 


+ - 


3  V 


*+1/2 


3  V 


*+1/2 


+ - 

3  V 


*+1/2 


e*+i/2  etc.  are  the  velocity  strains  in  the  calculation  of  the  stress  deviators. 
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Equivalent  Plastic  Strain 

{(«_  " p4 f  +  (p4  - pi. )’  +  («.  -  Pi.  f 
+  |[(P£.)!  +  (P£„)!+(P£„F]}  • 

Flow  Stress 

r+,=a[&  +  (£')“+1]\ 

Here  a,  b  and  c  are  material  constants,  not  to  be  confused  with  the  vector  components 
aijjc  etc. 


Time  Step 


(A/)"+3/2  =  0.67 


rn+l 


777? 


Min.  of  all  zones 


and  (Ar)"+3/2  ^l.lAr 


,*+1/2 


L  is  the  minimum  zone  thickness,  defined  as 


,*+1 


L"+1  =  -7^-,  where  vn+1  =  volume  of  zone  associated  with  point  i  jjc  at  tn+I,  and 


j*+1  is  the  area  of  the  largest  side  of  the  zone.  Also,  in  this  equation  for  At, 


a  =  sound  speed  calculated  from  the  equation-of-state  and 
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b=8[c2+cL]r+1 


•y  uo 

where  Q  and  Cl  are  the  quadratic  and  the  linear  q  constants,  respectively,  and  —  is  the 

dt 

rate  of  strain  used  in  the  calculation  of  q. 


Further,  (At)H+l  =  -(Ar‘+3/2  +  Ar"+,/2). 

Zt 


Boundary  Conditions 

Pseudo  zones  with  zero  mass  are  assumed  to  surround  the  grid  that  defines  the 
physical  object.  Thus  points  associated  with  the  surface  of  the  physical  object  may  be 
calculated  without  changing  the  logic.  Normally  a  free  surface  boundary  condition  is 
provided,  i.e.  the  pseudo  zone  pressures  are  considered  always  equal  to  zero.  Pressure 
boundary  conditions  may  be  applied  by  entering  the  desired  space-time  values  into  the 
pseudo  zones. 

A  reflection  boundary  condition  is  obtained  by  setting  equal  to  zero  the  normal 
component  of  accelerations  of  a  surface  point  when  it  points  into  the  reflection  surface. 


CHECK  PROBLEMS 
Simple  Harmonic  Motion 

The  calculation  of  the  motion  of  a  vibrating  plate,  clamped  at  one  end,  provides  a 
problem  that  can  be  readily  checked  by  elasticity  theory.  Orienting  the  plate  at  an 
arbitrary  angle  in  three  dimensional  space  activates  all  six  components  of  the  stress 
tensor. 


In  the  calculations  shown  in  Fig.  4  an  elastic  plate  clamped  at  the  top  is  set  into 
motion  by  applying  a  velocity  v  =  10  m/s  to  the  lower  right  edge  in  the  direction 
perpendicular  to  the  edge  for  a  time  t  =  50  ps.  After  this  time  the  applied  velocity  is 
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released,  but  the  lower  portion  of  the  plate  continues  to  move  due  to  the  kinetic  energy. 
Actually  upon  release  the  end  of  the  plate  initially  moves  faster  than  the  applied  velocity 
since  this  velocity  does  not  correspond  to  the  natural  frequency  of  the  plate.  Figure  4d  is 
a  time-displacement  plot  for  a  position  in  the  geometric  center  of  the  bottom  plane  of  the 
plate.  It  is  easily  verified  that  the  calculation  reproduces  the  fundamental  frequency  of 
the  plate. 

length :  L  =  52.5  mm 

Dimensions-  width :  W  =  20.0  mm 

thickness:  T  =10.0  mm 

bulk  modulus :  k  =  188  GPa 
Elastic  constants :  -  shear  modulus :  p  =  8 1 .4  GPa 

density :  p°  =  7.72  Mg/m3 

The  tensor  artificial  viscosity  used  in  this  calculation  is  Cns  -  0.05,  more  than  enough  to 
suppress  grid  oscillations  that  would  otherwise  occur.  Figure  4d  shows  that  the 
amplitude  of  the  oscillation  has  not  been  damped  or  affected  by  the  artificial  viscosity. 

Plasticity 

The  impact  of  a  right  circular  cylinder  on  a  rigid  boundary  provides  a  calculation 
to  test  the  plasticity  aspect  of  the  computer  program.  Since  this  problem  requires  only 
two  space  dimensions  it  can  be  calculated  with  the  HEMP  program.1  Figure  5a  shows 
results  of  the  HEMP  calculation  where  cylindrical  symmetry  is  incorporated  into  the 
fundamental  equations.  Figure  5b  shows  results  of  the  same  problem  calculated  with  the 
HEMP  3D  program  described  here.  It  can  be  seen  in  Figure  5b  that  the  cylinder  has  been 
discretized  with  three  dimensional  zones.  The  calculated  time  to  stop  the  cylinder,  30  ps, 
and  the  final  cylinder  length,  19.28  mm,  were  the  same  for  both  HEMP  and  HEMP  3D. 
Comparison  of  the  cylinder  profiles  at  t  =  30  ps  also  showed  almost  identical  results. 


27 


SLIDING  SURFACES  IN  THREE  DIMENSIONS 


The  sliding  surface  technique  described  here  has  evolved  over  several 
years  of  applications.  Very  good  results  are  obtained  even  for  severely  warped 
surfaces. 

The  implementation  of  sliding  surfaces  in  a  three  dimensional  Lagrange 
grid  i,  j,  k  follows  similar  procedures  as  slide  lines  in  the  two  dimensional 
problem.1  However,  instead  of  mapping  stresses  from  one  side  of  the  interface 
to  the  other  side  the  vector  accelerations  are  added  from  one  side  of  the  interface 
to  the  other.  (This  method  can  also  be  used  in  two  dimensions  but  there  is  no 
particular  advantage).  Interfaces  are  defined  in  i,  j,  k  space  that  separates  two 
regions.  The  grid  points  at  the  interface  of  one  region  slide  on  the  surface 
provided  by  the  grid  points  of  the  opposite  region  and  vice  versa.  The  grid 
points  associated  with  one  side  of  the  interface  are  designated  in  advance  as 
slave  points  while  the  grid  points  associated  with  the  opposite  side  of  the 
interface  are  designated  master  points.  The  calculations  are  symmetric  in  that 
the  grid  points  of  both  regions  at  the  interface  are  advanced  in  time  in  the  same 
manner.  After  the  grid  points  associated  with  each  region  have  been  advanced 
by  the  integration  time  step,  the  position  of  slave  points  are  adjusted  to  lie  on  the 
surface  defined  by  the  master  points  when  penetration  of  one  grid  surface  into 
the  opposite  grid  surface  occurs.  It  has  been  found  convenient  to  define  a  local 
surface  at  each  grid  point  as  the  plane  through  the  grid  point  that  is 
perpendicular  to  the  normal  vector  defined  at  the  point.  Thus,  the  interface 
between  two  regions  is  actually  composed  of  a  series  of  local  surfaces. 

All  grid  points  at  the  interface  of  the  two  regions  are  tagged  as  either  void 
open  or  void  closed.  Void  open  means  there  is  a  void  between  the  point  and  the 
opposite  surface  and  void  closed  means  the  point  is  in  contact  with  the  opposite 
surface.  Void  open  points  are  advanced  in  time  with  the  usual  free  surface 
calculations. 

At  the  interface  between  the  two  regions  it  will  be  convenient  to  refer  to 
the  point  that  is  currently  being  advanced  as  the  "current"  point.  Parameters 
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associated  with  the  other  side  of  the  interface  that  are  required  to  advance  a 
current  point  are  identified  by  the  word  "opposite."  The  roles  are  then  reversed 
after  calculations  have  been  completed  for  one  side  of  the  interface.  The 
symmetry  of  the  calculation  permits  sliding  surfaces  to  be  defined 
simultaneously  in  more  than  one  direction.  However,  for  illustration  of  the 
method,  we  will  assume  a  sliding  surface  at  a  constant  Lagrange  coordinate  j;  see 
Fig.  6.  The  letter  f  will  be  used  to  designate  a  current  point. 

Free  surface  boundary  conditions  are  used  to  calculate  the  acceleration  of 
point  f  in  the  x,  y,  z  coordinate  system.  The  components  of  acceleration  are 
transformed  into  a  coordinate  system  where  two  components  are  in  die  plane  of 
the  sliding  surface  interface  at  point  f.  The  acceleration  component  normal  to  the 
interface  includes  a  contribution  of  mass  from  the  opposite  grid  and  in  addition, 
the  normal  component  of  acceleration  of  the  opposite  grid.  The  two  acceleration 
components  in  the  plane  of  the  interface  are  unchanged  by  the  presence  of  the 
interface.  The  normal  component  of  acceleration  from  the  opposite  grid  must 
include  a  contribution  of  mass  from  the  present  grid.  Thus,  the  symmetric 
treatment  of  the  interface  calculations  requires  preprocessing  each  side  of  the 
interface.  A  final  calculation  is  then  made  to  advance  in  time  points  associated 
with  each  side  of  the  interface. 


OUTLINE  OF  SLIDING  SURFACE  CALCULATION 


Step  I 


For  every  grid  point  on  the  sliding  surface  interface,  calculate  in  advance 
the  following  quantities  and  store  with  the  grid  point. 

1 .  The  mass  per  unit  area,  m. 

2.  The  acceleration  assuming  the  point  is  on  a  free  surface. 

3.  An  outward  pointing  unit  vector  normal  to  an  element  of  surface 
defined  at  each  point  on  the  sliding  surface  interface. 
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Step  n 


For  each  current  void  closed  grid  point  f: 

1.  Locate  the  three  points  of  the  opposite  grid  that  are  closest  to  point 
f. 

2.  Determine  by  interpolation  the  mass  per  unit  area  of  the  opposite 
grid  at  the  position  of  point  f. 

3.  Calculate  the  mass  weighting  factor  at  point  f  (z-factor). 

4.  Resolve  the  acceleration  obtained  in  Step  1-2  in  the  direction  of  the 
normal  vector  defined  in  Step  1-3. 

5.  Repeat  1-4  with  the  opposite  grid  as  the  current  grid. 

Stepm 

For  each  current  void  closed  grid  point  f: 

1.  Locate  the  three  points  of  the  opposite  grid  that  are  closest  to  point 
f. 

2.  Determine  by  interpolation  a  value  for  the  normal  component  of 
acceleration  from  the  opposite  grid  at  the  position  of  current  grid 
point  f. 

3.  Add  to  the  normal  component  of  acceleration  of  the  present  grid 
from  step  n-4  the  normal  component  of  acceleration  contributed  by 
the  opposite  grid  from  Step  2  above. 

4.  Calculate  new  velocities  and  new  coordinates  in  the  x,  y,  z 
coordinate  system. 
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5.  Repeat  1-4  with  the  opposite  grid  as  the  current  grid. 

Step  IV 

1.  Test  for  penetration  of  the  current  grid  points  into  the  opposite 
grid. 

2.  Adjust  the  velocities  for  all  points  that  have  penetrated  the  opposite 
grid  to  conserve  momentum. 

3.  Repeat  (1)  and  (2)  with  the  opposite  grid  as  the  current  grid. 

4.  Relocate  slave  grid  points  that  have  penetrated  the  master  grid  onto 
the  master  grid.  Void  open  points  are  advanced  in  time  in  the  x,y,z 
coordinate  system  using  free  surface  boundary  conditions 
independent  of  Steps  I,  II  and  HI.  If  penetration  of  the  opposite  grid 
occurs  from  Step  IV  the  point  is  relocated  and  labeled  void  closed. 
Void  closed  points  are  labeled  void  open  when  the  distance  of  the 
point  to  the  opposite  surface  is  greater  than  10%  of  the  zone  size. 


Calculational  Steps  to  Advance  in  Time  Grid  Points  on  a  Sliding  Surface 
Step  I 

1.  Calculate  the  mass  per  unit  area  for  all  grid  points  on  the  sliding 
interface.  Referring  to  Fig.  6  assume  point  a  is  a  point  on  the 
interface.  The  mass  per  unit  area,  ma,  is  given  by: 

A/0  +  M9  +  M9  +  M9 
Ad  +  As  +  As  +  A» 


M©  is  the  mass  of  zone  ©  etc.  A©  is  the  area  of  the  triangle  in  zone  © 
that  is  associated  with  opposite  grid  point  a.  Similar  for  A®  ,  A®  ,  and 
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A@  .  The  parameter  ma  is  seen  to  be  the  average  of  the  mass  per  unit  areas 
of  the  zones  that  share  point  a. 


2.  Calculate  the  acceleration  of  each  point  on  the  interface  with  free  surface 
boundary  conditions. 

a.  For  a  given  point  i,  j,  k  calculate  the  acceleration  A*j  k  (see  point  a 
Fig.  6). 


dx  r 

—  i 

dt 


dy  ?  dz 

+  —  j  +  — 

dt  dt 


k 


i)  x  direction 


Pi.J.k  &  *  J  i.JJt 


,  where 


{(^«)®[(}V/  3vXziK  ~  zv)  (zvi  ~  zv)(y/v  —  >v)] 


+  (^**)®[(y//  yvtzvi  zv)  (zn  zv){yvl  jv)] 

+  -y,„){zvi  ~  ziii)~{ziv  -z„,){yvi  ~ yin )] 

+  (^«)o[(^w  y  m  )(z//  ~  zm)~  (zv/  “  zm  X^// ~  ?///)]} 

(4>)y  ..=|K+W*  +  M.+M.] 

To  form  (l/p  ctT^/dy)*^  ,  replace  each  I,**  in  the  right  side  of  the  above 

expression  with  Txy,  every  y  with  the  corresponding  z,  and  each  z  with  the 
corresponding  x. 

To  form  (l/p  cTu /dz)*. k  ,  replace  each  ZxX  in  the  above  expression  with 
Tzx,  every  y  with  the  corresponding  x,  and  every  z  with  the  corresponding  y. 
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ii)  v  direction 


(ur„x  fiaz„Y 

— r2-  =  same  as H2- 

U  ^  Am  Ip  *  Am 


defined  above. 


except  replace  each  Xxx  by  the  corresponding  value  of  Txy. 


fiax.Y  ri»;Y 

- r-2-  =  same  as  — ~ 

VP  %  A.m  Ip  ❖  Am 


defined  above. 


except  replace  each  Txy  by  the  corresponding  value  of  lyy . 


fiar*Y  fiar„Y 

— -2.  =  same  as  — r2- 

\P  *  kit  U  *  kit 


defined  above. 


except  replace  each  Tzx  by  the  corresponding  value  of  Tyz 

iii)  z  direction 


{* Y  -  1  fa  fg* 

Ur  Ay.*  P"y.t  L<*  ^ 


<?Zn 
+  A  . 


=  same  as 


ip  *  A,,  tp*. 

except  replace  each  Lxx  by  the  corresponding  value  of  Tzx. 


jggfcY 

P  *  J,,t 


defined  above. 


1  ^  \* 

— ~  =  same  as  — r2-  defined  above, 

p  ❖  A/..  Ip  ❖  A/.. 


except  replace  each  Txy  by  the  corresponding  value  of  Tyz. 
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defined  above, 

i.j.k 

except  replace  each  Tzx  by  the  corresponding  value  of  Zzz- 

Determine  an  outward  pointing  unit  vector  normal  to  an  element  of 
surface  defined  at  each  point  on  the  sliding  surface  interface.  Referring  to  Fig.  6 
assume  point  a  is  any  point  on  the  sliding  surface  interface. 

a)  Calculate  the  normal  vectors  for  each  of  the  triangular  surface  areas 
associated  with  point  a.  The  normal  vector  corresponding  to  zone 
©is: 


.P  * 


=  same  as 


i#H 

p  dz 


i 


r*y.tv  ~ 


xv-xa 
xiv  ~xu 


j 

yv-ya 

y/v-y. 


~  ^a.V.lV  1  +  Ba.VJV  J  +  Cay. IV  k 


where 


\y.N  =  [(y^ — y«)  (z/v  -  z«) — {zv  ~  z«)  {y/v — y«)] 
y,IV  =  *[(^V  ~Xa)  (Z/V  —  Z«)~  {ZV  ~  Zm)  {XIV  ~  •*«)] 

Ca,v,iv  =[(*v  ~  x a)  (y/v  ~ y«)— (yv —  y«)  C*/v — •*«)] 


Similar  for  zones  ®,  ®,  and©.  Fig.  6. 

Note:  The  vector  cross  products  must  be  taken  so  that  the  normal 
vectors  point  outward  from  the  grid.  A  convenient  way  to  assure 
that  vector  ra  VJV  points  outward  from  the  grid  is  to  take  the  dot 
product  of  raVlv  with  the  vector  formed  by  point  a  and  its 
corresponding  interior  point.  The  interior  point  is  called  a 
connector  point.  In  Fig.  6  point  VI  is  the  connector  point  for  point 
a.  If  the  product  is  positive  reverse  the  sign  raVJV  otherwise  the 
vector  has  the  correct  outward  direction. 
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b)  Calculate  ra,  a  unit  vector  obtained  from  the  average  of  the  vectors 
perpendicular  to  the  triangles  that  surround  point  a,  i.e.,  Fig.  6 
triangles  (a,V,IV),  (a,IV,in),  (a,!!!,!!)  and  (a,n,V). 


ra  =  AJ  +  Ba]  +  Cjc  /^A2  +  B2  +  C 2  =li  +mj  +nk 


where  Aa/  Ba/  and  Ca  are  the  direction  numbers  of  vector  ra  and  1, 
m  and  n  are  the  direction  cosines. 


\  —  (Ai.v./v  +  Aajv,in  +  \,m.u  "*■  A a.II.V ) 

Ba  =  (BaVrv  +  Ba,IV.III  +  Ba 

MUi  +  Ba.ll.v) 

Ca  =  {Ca.V.lV  +  -'ll. IV, III  +  Ca.1,1.11  +  Cajiy  ) 


StepII 

1.  Locate  the  opposite  surface  points  associated  with  each  current  void 

closed  point  f: 

a)  Calculate  d2,  the  square  of  the  distance  from  point  f  to  successive 
points  i,  k  of  the  opposite  grid. 

d) = W -■ *<.*  f + Ov  -  yJf + ( v -  v )’ 

b)  Let  point  a,  of  the  opposite  grid  be  the  point  which  has  the  shortest 
distance  from  point  f.  See  Fig.  7. 

c)  Project  the  points  1,  2,  3,  4  and  f  onto  the  local  surface  at  point  a, 
defined  by  point  a  and  the  unit  vector  ra  at  point  a.  Designate  as 
(x*  y*,  z*)i  the  coordinates  of  a  point  i  and  (x,  y,  z)j  the  coordinates 
after  projection  onto  the  local  surface  at  point  a. 


x,  =  x-  -  ldi 
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y, =y‘-mdt 

z, .  =  z*  -  ndi 

where  di=ra  •  ci  1  =  1,2, 3, 4,  and  f.  The  vectors  ci  are  formed  by  the 
points  i  and  the  origin.  Here  /,  m  and  n  are  the  direction  cosines  of 
ra,  the  unit  vector  at  point  a  calculated  in  Step  I.3b. 

d)  Point  f  can  be  in  any  one  of  the  four  quadrants  formed  by  the 
projections  of  opposite  grid  points  surrounding  point  a  onto  the 
surface  defined  by  point  a  and  the  unit  vector  at  point  a.  Fig.  7.  The 
following  procedure  is  used  to  locate  the  quadrant  that  contains 
current  point  f . 

i)  Referring  to  Fig.  7  calculate  the  area  of  triangles  Aa<2,3,  a,  p, 
and  y  using  the  coordinates  obtained  in  Step  lc  above. 

ii)  Point  f  is  contained  within  triangle  Ag^  if:  [Aa/2,3  -  (a  +  P  + 
y)  <  Kb3  Aa  2,3].  Here,  Ag^,  a,  P,  y  refer  to  the  areas  of  the 
triangles. 

in)  Repeat  step  (2)  for  the  remaining  quadrants. 

iv)  If  all  quadrant  tests  fail  to  locate  point  f  then  select  the  next 
closest  point  a  to  point  f  and  repeat  (1-3). 

v)  Repeat  (1-4)  for  all  points  f  in  the  current  plane. 

Note:  It  is  important  that  the  search  logic  described  above  be 
conducted  with  the  grid  points  projected  onto  the  local  surface 
(Step  lc). 

e)  It  is  necessary  to  determine  when  a  current  point  f  is  not  covered  by 
the  opposite  grid  as  shown  in  Fig.  8.  Assume  current  point  f  has 
been  determined  closest  to  opposite  grid  point  a  and  it  is  known 
that  grid  point  a  is  on  the  boundary  of  the  opposite  grid.  Fig.  8.  The 
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four  quadrants  for  the  search  routine  described  above  are 
formulated  by  extending  the  opposite  grid  through  point  a. 
Referring  to  Fig.  8  an  extension  point  3  is  established  by  calculating 
a  vector  extension. 

i)  v,  =  v.  +  k.[v.-v,) 

Here  Va  designates  the  vector  formed  by  point  a  and  the  origin. 
Similar  for  points  1  and  3.  The  parameter  ke  provides  the 
dimension  of  the  extension.  kt  is  taken  as  a  large  number,  e.g., 
1000,  to  assure  that  point  f  is  covered  by  the  opposite  grid  for  the 
search  routine  that  locates  the  quadrant. 

ii)  Coordinates  of  extension  point  3  (Fig.  9) 

xi=x,+k'(xa-xl) 

yi=y.+k.(ya~yi) 

Z3  =  Z«+*«(Z.-2l) 


iii)  When  point  f  has  been  located  in  a  quadrant  that  contains 
the  extension  point  it  is  outside  the  opposite  grid.  If  the 
point  is  more  than  one  zone  thickness  off  the  opposite  grid  it 
is  considered  a  free  point  independent  of  the  sliding 
interface.  If  the  point  is  less  than  a  zone  thickness  off  the 
opposite  grid  it  is  considered  still  on  the  sliding  interface.  To 
make  this  distinction  the  extension  point  is  recalculated  with 
a  value  of  kt  to  provide  an  extension  of  the  opposite  grid  of 
approximately  one  zone  thickness  of  the  current  grid.  (The 
default  value  is  kt  =  1,  which  assumes  both  grids  are  the 
same  size). 

Point  e  in  Fig.  8  is  the  new  extension  point.  The  following 
method  is  used  to  locate  the  position  of  current  point  f  with 
respect  to  the  extension  surface.  Calculate: 
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Here  Ae,  Be  and  Ce  are  the  direction  numbers  of  vector Cr . 
S  is  the  vector  from  point  a  to  current  point  f  and  Ct  the 
vector  from  point  a  to  the  extension  point  e. 

If  de  >  I,  point  f  is  beyond  the  extension  and  is  accelerated 
with  free  surface  boundary  conditions. 

If  de  <  I,  point  f  is  considered  to  be  still  on  the  sliding 
interface  and  the  calculation  proceeds. 

2.  Calculate  an  interpolated  mass  per  unit  area,  mf,  at  the  position 
corresponding  to  point  f.  Figure  9  shows  an  overlay  of  the  current  grid 
containing  point  f  on  the  opposite  grid.  We  wish  to  obtain  the  mass  per 
unit  area  of  the  opposite  grid  at  the  position  of  current  point  f.  This  mass 
per  unit  area  will  then  be  used  to  increase  the  mass  associated  with  point  f 
for  the  acceleration  of  point  f  in  the  direction  normal  to  the  sliding 
interface. 


Referring  to  Fig.  10  the  mass  per  unit  area  at  point  f  is: 

(m  )-  maa  +  mb(5+mcy 
'  f'  a  +  P+y 


See  Step  1-1  for  the  calculation  of  ma,  etc. 

3.  Calculate  the  mass  weighting  factor  at  point  f.  (z  factor) 

a)  Let  Mm  be  the  mass  due  to  the  opposite  surface  that  is  to  be 
included  with  the  mass  of  point  f. 

Mm  =  (mf)  [A®  +  A®  +  A®  +  A®]  (see  Fig.  9) 
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Here  A®  is  the  scalar  aiea  of  the  triangle  (f,  V,  IV)  in  zone  (D  Fig.  11, 
and  is  calculated  as  follows: 


~  ^(A/.v./v)  +(^/.v./v)  +(^7.v'./v) 


where 


A/.v.tv  -  [(jv  ~yj)  { zv  -  zf )  -  ( zv  ~z/)  (ytv  -  >7)] 

Bf.v.iv  —  [('XV'  ~  Xf)  (Z/V  —  Zf)  ~  [ZV  ~  Zf)  (XIV  ~  */)] 

Cf.v.iv  —  [(*v —  xf)  {yrv  ~y/)~{yv  ~ 37)  -■*/)] 

Similar  for  2A®  ,2A®and2A®. 
b)  Calculate  the  z-factor 


A  “  Jt  T-j - 

— (M&  +  Af«  +  +  M9) 

Here  M®  etc.  are  the  masses  associated  with  point  f,  see  Fig.  11. 

4.  Calculate  the  acceleration  of  grid  points  on  the  sliding  surface.  The 
acceleration  normal  to  the  sliding  surface  includes  the  mass  from  the 
opposite  grid  using  the  z-factor  determined  from  the  preceding  step. 

a)  Calculate  N/t  the  free  surface  acceleration  of  point  f  resolved  in  the 
direction  of  the  average  normal,  rf . 


The  free  surface  acceleration  of  point  f,  A 'f,  was  calculated  in  Step  I- 
2.  The  average  normal,  rf,  was  calculated  in  Step  I.3.b. 
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b)  Calculate  Nf,  the  acceleration  of  current  point  f  in  the  direction  of 
the  average  normal  that  includes  the  mass  of  the  opposite  grid. 


5.  Repeat  1-4  with  the  opposite  grid  as  the  current  grid. 

Note:  Nf  is  a  partial  acceleration  normal  to  the  surface  defined  at  point  f 
that  includes  the  mass  of  the  opposite  grid.  The  total  normal  acceleration 
of  current  point  f  must  also  include  a  contribution  from  the  opposite  grid 
and  is  described  in  Step  III  that  follows. 


Step  in 

1.  For  the  current  void  closed  point  f,  locate  the  opposite  grid  points  that 
surround  point  f.  See  Step  13-1  and  Fig.  9. 

2.  From  the  three  opposite  grid  points  a,  b,  c  that  surround  point  f.  Fig.  10, 
determine  acceleration  vector  N™.  Vector  N™  is  an  interpolated  normal 
component  of  acceleration  from  the  opposite  grid  at  the  position  of  point  f. 
The  interpolation  method  is  shown  in  Fig.  10. 

fjoc  ^Naa  +  Nbp  +  Nj 

1  a+P+Y 

Na,  Nb  and  Ne  are  the  accelerations  of  opposite  grid  points  a,  b  and  c 
calculated  in  Step  II-5. 

Test  for  void  opening.  If  Nf-Cf>-  0  and  Cf< 0  tag  point  f  as  void 
open.  Nf  is  calculated  in  Step  II-4.  The  vector  C{  is  formed  by  point  f  and 
the  connector  point  associated  with  point  f.  If  the  above  test  is  positive  the 
acceleration  of  point  f  is  given  by  A‘f,  the  free  surface  acceleration  from 
Step  1-2. 
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3.  Calculate  the  total  acceleration,  Af,  of  point  f. 

Af  =  A)  -  N}  +  Nf  +  N™  =  Aj  +  Aj  +  Ajc 

Note:  Nf  must  be  saved  for  use  in  the  interpolation  procedure  when  the 
above  process  is  reversed  and  the  current  grid  becomes  the  opposite  grid. 

4.  Calculate  the  x,  y,  z  components  of  velocity  and  new  coordinates  for 
current  point  f. 

a)  Velocity 

x}+U2=x}-in+AtmAx 
y,*in  =  y}~u2  +  At*  Ay 
z*/u2  =  i*f-U2+At*Ax 

b)  New  coordinates 

xf'=x*/+x}+U2At*+U2 
y*+l  =  y*  +yHf+U2A t**U2 

z;+,  =  z;+i;+,,2Ar"+1/2 

5.  Repeat  1-4  for  all  interface  grid  points. 


Step  IV 

1.  Test  to  see  if  a  point  f  has  penetrated  the  opposite  grid.  Assume  f,  in  Fig. 
7,  is  a  point  on  a  grid  that  is  to  be  tested  for  penetration  into  the  opposite 
grid  local  su!r?i  c  at  point  a. 

a.  Calcuwfc  d  the  perpendicular  distance  from  the  point  f  to  the  local 
surface  at  point  a. 
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d  =  [ii+m]  +  nk][(xf-xa)i+(yf-ya)]  +  (zf-za)k] 


or 


d  =  +[/(*/  ~x°)  +  m(yf  ~ y*)+  n(zf  ~  z-)] 

Here  /,  m  n  are  the  direction  cosines  of  the  unit  vector  defined  at 
point  a  (see  Step  L3b). 

b.  If  0  <  d  <  +  6,  point  f°+1  remains  as  calculated  in  Step  HI  with  the 
same  void  status  as  before.  Here  8  is  a  positive  number  equal  to  0.1 
times  the  grid  spacing  calculated  as: 

8  =  0A^A2  +  B2  +  C2 

c.  If  d  >  5,  point  f*1*1  remains  as  calculated  in  Step  HI  and  is  tagged 
void  open. 

d.  If  d  <  0  point  fn+1  has  penetrated  the  opposite  grid  and  is  tagged 
void  closed. 

2.  Adjust  the  velocities  of  all  void  closed  points. 

a.  Calculate  velocities  normal  to  the  interface.  Assume  point  f  in  Fig. 
9  has  penetrated  the  local  surface  at  point  a  of  the  opposite  grid, 
calculate  the  velocity  components  Nf,  Nc,  Nb  and  Ne  that  are 
normal  to  the  surface. 

Na  =  lxa+mya  +  nia 
Nb  =  lxb  +  myb  +  nzb 
Ne=lxc+  myc  +  nzc 
Nf  =  lx f  +  myf  +  nij 
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b.  Calculate  N+f  the  velocity  of  point  f  from  the  conservation  of  linear 
momentum: 


a  +  p  +  y  '  ' 

aM.  +  pMt  +  *i,  t 
a+p  +  y  1 


Mf  +  j[M9  +  M.  +  M*  +  Mm]  =  2(0), 

Ma  +  +  M9  +  M9  +  Ma]  =  2(0). 

similar  for  Mb  and  Mc. 

Note:  The  mass  C>  associated  with  a  point  on  the  interface  is 
calculated  in  Step.  L 

Note:  To  minimize  the  number  of  calculations,  only  the  linear 
momentum  has  been  considered  instead  of  including  conservation 
of  angular  and  linear  momentum  as  was  done  in  the  two 
dimensional  problem.  Actually  it  is  the  artificial  viscosity,  q,  and 
the  equations  of  motion  that  accomplish  the  conservation  of 
momentum.  Adjusting  the  velocities  at  the  interface  after  a 
collision  sets  up  the  initial  conditions  for  the  artificial  viscosities  on 
each  side  of  the  interface. 

c)  Calculate  the  x,  y,  z  components  of  the  new  velocity  of  point  f. 

x?U2=xf+lAN, 
yfU2=yf+mANf 
zH/u2  =  zf  +  nANf 


where 
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a*/ = (*;-*/) 

Here  x),  zf  refers  to  the  velocity  of  point  f  before  the  adjustment 
for  conservation  of  momentum.  The  velocities  of  all  void  closed 
points  are  adjusted  point  by  point  That  is,  only  the  velocity  of  the 
point  under  consideration  is  adjusted,  point  f  in  the  example 
described  above.  The  velocities  of  all  mass  points  on  one  side  of  the 
grid  are  adjusted  when  penetration  has  occurred.  Subsequently  the 
velocities  of  all  points  on  the  opposite  grid  are  adjusted.  Thus,  all 
of  the  old  velocities  and  coordinates  must  be  retained  until  all  of 
the  velocities  of  both  sets  of  sliding  surface  points  have  been 
adjusted. 

3.  Declare  one  grid  the  slave  grid  and  the  grid  opposite  it  the  master  grid. 
Relocate  slave  points  onto  the  master  surface  for  slave  points  where  d  <  -6. 
A  slave  point  f  that  has  penetrated  the  master  surface  is  set  back  to  the 
master  surface  by  subtracting  the  length  d  from  the  position  of  the  point. 
For  the  sign  convention  used  here  d  is  a  negative  number.  The  direction 
cosines  /,  m,  n  point  outward  from  the  grid,  thus  the  new  coordinates  of 
point  f  are: 

x}+l=x'f-ld 

x/+1  -y/~m d 

z?'  =  z'f-nd 

Here  xf,  y‘r  z‘f  refers  to  the  coordinates  of  point  f  used  to  determine  that 
d<-8. 


When  penetration  occurs  new  velocities  are  calculated  on  both  sides  of  the 
interface,  but  only  the  position  of  slave  points  are  adjusted. 
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APPLICATIONS  OF  SLIDING  SURFACE  ROUTINE 

The  major  difficulty  with  sliding  surface  routines  arises  from  failure  of  the 
search  routines  that  must  locate  one  grid  with  respect  to  the  other.  The  problem 
becomes  aggravated  with  curved  or  warped  surfaces  with  a  search  routine  that 
operates  in  three  dimensions.  The  method  used  here  projects  the  grid  point  onto 
local  two  dimensional  surfaces  to  establish  the  orientation  of  one  grid  with 
respect  to  the  other.  With  this  procedure  the  search  technique  is  robust  even  for 
distorted  surfaces.  Figure  12  shows  an  application  with  two  curved  surfaces. 
Fig.  13  shows  the  acceleration  of  a  metal  plate  by  an  explosive  with  a  sliding 
surface  between  the  two  materials.  The  explosive  was  detonated  at  nine  equally 
spaced  points  on  a  line  along  the  top  surface  of  the  explosive.4 

ZONE  DIMENSION  CHANGE  AND  SUBCYCLING 

It  is  useful  to  be  able  to  change  from  coarse  to  fine  zoning  in  a  localized 
region  and  to  be  able  to  join  two  independent  grids  at  an  interface.  The  latter 
being  especially  important  for  constructing  grids  for  three  dimensional  problems. 

Zone  dimension  change  at  an  interface  in  two  dimensions 

Figure  14  shows  schematically  a  zone  change  from  large  to  small  zones 
across  a  Lagrange  coordinate  kt.  The  grid  with  the  largest  zone  size  is  chosen  as 
the  master  grid  and  the  small  zone  grid  the  slave  grid.  The  master  grid  defines 
the  interface  k,.  In  Fig.  14  grid  points  associated  with  the  master  and  slave  grids 
are  shown  as  closed  circles  and  open  circles,  respectively. 

1.  To  advance  in  time  master  and  slave  points  on  the  interface  kt,  Fig.  14. 

a.  Calculate  the  partial  acceleration  of  all  master  points  j,  k  on  k-line 
k,. 


Components  of  partial  acceleration  for  point  a.  Fig.  14. 
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Oa  =  " y'"} + -  y») 

~{T*y  )®  (x«  ~  *7// )  “  {T*y  )®  (xff/  ~  xPf  )] +  ^ "  (a)/.*  ]• 

~[Tv )®  (?«  ~  y«i )  ~  {T*> )®  {y*u  ~  y* )] + ^**(0)^ } 


The  z  factor  that  is  found  in  equations  (i)  and  (ii)  above  is  obtained 
by  mapping  the  masses  of  the  slave  zones  between  master  points  II 
and  IV  onto  the  master  grid. 


b.  Calculate  the  partial  accelerations  for  all  slave  points  on  k-line,  kt. 
The  same  procedure  as  above  is  used.  The  z  factor  now  maps  mass 
from  the  master  grid  onto  the  slave  grid. 


c.  For  each  master  grid  point  on  k-line  k ,  determine  a  slave  grid 
partial  acceleration  by  interpolation.  Referring  to  Fig.  14  the  slave 
grid  partial  acceleration  corresponding  to  master  grid  point  a  is: 


+  (l-a) 


Here  a  =  lat/l„  where  l„  is  the  distance  between  points  a  and  t  and 
/„  the  distance  between  points  s  and  t. 


Similar  for 
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d.  Calculate  the  total  acceleration  of  all  master  points  on  k-line  kg. 
Referring  to  Fig.  14  the  total  acceleration  of  master  point  a  (point  j, 
k  Fig.  14)  is: 

l dt)j +  { dt)a  +{dt)a 

Similar  for 


e.  Calculate  new  velocities  for  all  master  points  on  k-line  kg. 


f.  Obtain  new  velocities  for  the  slave  grid  points  on  k-line  kg  by 
interpolation  so  that  the  original  spacing  between  consecutive 
master  points  is  maintained.  Referring  to  Fig.  14  the  new  velocities 
for  slave  point  t  are: 

P)x?U2 
P)y?in 


y?'12  =  py?m  +{ i- 
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Here  p  is  a  constant  calculated  when  the  grid  is  generated.  4,  is  the 
distance  from  point  t  to  point  b  and  the  distance  between  points 
a  and  b. 


Zone  dimension  change  at  an  interface  in  three  dimensions 

The  same  procedure  is  followed  as  described  for  the  two  dimensional 
case.  Fig.  15  shows  two  grids  that  are  to  be  joined  together  without  the  interface 
grid  points  of  both  grids  being  necessarily  coincident. 

1.  To  advance  in  time  master  and  slave  points  at  the  interface  of  two 
grids. 


a.  Calculate  the  partial  acceleration  for  all  interface  grid  points 
associated  with  the  master  grid.  Referring  to  Fig.  15  we  wish 
to  accelerate  grid  point  a  associated  with  master  grid  zones 
©,  ®,  ®,  ®. 

Partial  acceleration  in  x-direction. 


i) 

ii) 


fdxX^  1 
\dtJij*  Pi,J,k  . 


dx  dy 


arT 

dz  .  ’ 


where 


p(  <&*) .  ,k  ~  z4d>.  -k  yv )  ( Zlv  *v)-(*v,-*vblv  ~>v)] 


+  (£«  )®[(^//— ^v)  (zv,  -  zv)~{zu  -  zv)bv,  —  3V)] 

+  -yin)  {zvi  ~zm)~{ziv  —  zin){yvi  ~ ym)\ 

+  (^»)o[0V/  —  7///)  (Z1I  ~  zm )  ~  (ZVI  —  ZUI )(>'//  -  3Vff)] 

Mi)  (%  =  i[M„  +  M®  +  Ma  +  Me] 
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The  factor  z  that  appears  in  (ii)  above  is  the  weighting  factor  that  maps  the 
mass  of  the  opposite  grid.  The  remaining  terms  in  equation  (i)  are  composed  in 
the  usual  manner  and  include  the  z  factor,  as  shown  in  equation  (ii),  the  finite 
difference  equation  for  the  first  term  in  equation  1.  In  a  similar  manner  the  y,  z 
components  are  calculated. 

b.  Calculate  the  partial  acceleration  for  all  points  associated 
with  the  slave  grid  on  the  interface.  The  z-factor  now  maps 
the  mass  from  the  master  surface  onto  the  slave  grid-point. 

c.  For  each  master  grid  point  determine  a  slave  grid  partial 
acceleration  by  interpolation.  Referring  to  Fig.  16  assume 
master  grid  point  a  has  been  found  to  be  in  the 
neighborhood  of  slave  grid  points  f,  g,  h. 

f.r  -(fHSHS). 

\dt)m  a+p+y 

Here  f,  g,  and  h  denotes  the  partial  acceleration  of  slave  grid 
points.  Similar  for  the  y,  z  components. 

d.  Calculate  the  total  acceleration  for  all  master  points  on  the 
interface.  For  master  point  a  (point  i,  j,  k.  Fig.  15)  the  total 
acceleration  is: 

(dn  f  #  Y1™  f  di  Y“,“r 

V  dt  Ji.y.t  l  dt)  +  l  dt) 

Similar  for  y,  z  components  of  acceleration. 

e.  Calculate  new  velocities  for  all  master  points  on  the 
interface. 

f.  Obtain  new  velocities  for  the  slave  grid  points  on  the 
interface  using  the  interpolation  scheme  above. 
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Subcycling  with  zone  dimension  change  in  two  dimensions 

The  time  step  for  a  given  cycle  is  dictated  by  the  zone  with  the  smallest 
zone  dimension  divided  by  the  local  sound  speed.  Rather  than  calculate  the 
entire  grid  with  a  given  time  step  a  saving  in  computer  time  can  be  obtained  by 
dividing  the  grid  into  different  regions  with  a  different  time  step  for  each  region. 
A  region  with  small  zones  is  calculated  for  several  time  steps  until  the  time  step 
of  a  region  that  can  use  a  larger  time  step  is  reached.  The  region  with  the  larger 
time  step  is  then  advanced  with  a  single  time  step  equal  to  the  sum  of  the  time 
steps  used  in  the  region  with  the  small  zone.  It  is  convenient  to  calculate  the 
region  with  the  largest  time  step  first.  The  time  steps  for  the  regions  with  the 
smaller  time  step  requirements  are  chosen  so  that  an  integral  number  of  equal 
time  steps  can  be  used  to  reach  the  time  step  used  for  the  largest  grid. 


Example  for  a  zone  size  change  of  two  to  one 

1.  With  a  time  step  At  calculate  new  velocities  for  all  grid  points  of  the 
largest  grid  using  the  stress  boundary  conditions  provided  by  the 
small  grid.  Advance  all  points  of  the  large  grid  and  calculate  the 
new  zonal  parameters.  Save  the  old  positions  of  the  large  grid  for 
points  on  the  interface  between  the  two  grids. 

2.  Find  the  new  velocities  of  all  small  grid  points  that  are  on  the 
interface  by  interpolation.  These  velocities  are  the  boundary 
conditions  during  the  subcycling  of  the  small  grid. 

3.  With  the  old  interface  positions,  so  that  all  of  the  positions  of  the 
small  grid  points  are  at  the  same  time  calculate  new  velocities  from 
the  acceleration  equations  for  all  small  grid  points  except  those  on 
the  interface  using  a  time  step  At/2. 

4.  Calculate  new  coordinates  for  all  of  the  small  grid  points  including 
those  on  the  interface.  The  points  on  the  interface  use  velocities 
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from  step  b  to  obtain  new  coordinates  and  the  remainder  of  the 
points  of  the  small  grid  used  in  the  velocities  from  Step  c. 

5.  Calculate  the  new  zonal  quantities  for  the  small  grid. 

6.  Calculate  new  velocities  again  for  all  small  grid  points  except  those 
on  the  interface  using  time  Step  At/2. 

7.  Calculate  new  coordinates  for  all  grid  points  including  those  on  the 
interface.  (The  large  grid  coordinates  will  now  coincide  with  the 
values  obtained  in  Step  1. 

8.  Calculate  new  zonal  quantities  for  the  small  grid. 

The  procedure  is  similar  for  three  dimensions. 
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(i,  j+l/k+l) 


(i+  1,  j+  1,  k+  1) 


(•+  1/  j+  1#  k) 


Fig. 


Grid  numbering  scheme  for  Zone  ©. 


Fig.  2. 
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*  1,  I  1  11  I  1.1  |  I  1  -  1  —A  1 


Simulation  of  the  motion  of  a  vibrating  elastic  plate,  (a)  Position  of 
maximum  positive  displacement,  t  =  275  ns.  (b)  Position  of  maximum 
kinetic  energy,  t  =  360  (is.  (c)  Position  of  maximum  negative  displacement, 
t  =  450  ps.  (d)  Displacement  history  for  a  point  in  the  geometric  center  of 
the  bottom  plane. 


Fig.  4. 
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a  current  point  f. 

(a)  Quadrants  surrounding  opposite  grid  point  a  with  current  grid  point  f 
in  quadrant  3. 

(b)  System  of  triangle  shown  for  quadrant  3. 


Fig.8.  Schematic  grid  to  determine  if  a  current  point  f  is  outside  the  opposite 
grid. 


58 


current  grid - 

opposite  grid - 


Fig.  9.  Current  point  f  is  associated  with  the  opposite  surface  formed  at  point  a. 


a  =  area  of  Acbf 
(3  =  area  of  Acfa 

y=areaof  Afba 


m/  = 


a+f)  +  y 


ma,  m^,  mc  are  the  mass  per  unit 
area  of  opposite  grid  points  a,b,c. 


Fig.  10.  Weighting  scheme  for  obtaining  the  value  of  a  parameter  defined  at 
points  a,  b,  c  at  position  f. 
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Steel  case 


High  explosive 


Copper  liner 
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ig.  13.  Calculation  of  a  three-dimensional  implosion  of  a  copper  liner. 


(a)  Section  view  of  geometry. 

(b)  Time  sequence  of  implosion. 


Fig.  14.  Schematic  of  a  zone  dimension  change  at  k-line  ks. 


Fig.  16.  Interpolation  scheme  for  obtaining  at  position  a  information  defined 
positions  f,  g,  h. 


